For a multi-site MRI study on MS, we want to understand the effect of a clinical measure, EDSS, on imaging metrics. Generally the distribution of EDSS is not normal, but is bimodal normal, where scores cluster around 2(low severity) and 6(high severity) on a scale of 0 to 10. The goal of this notebook is to first group subjects into low and high severity groups based on EDSS scores. We will then run the simulated data through a Logistic Regression to predict group membership.
In order to get the most power behind the coefficient estimate for EDSS, we need to combine data across the multiple sites. When combining data, we aim to calibrate the scanner metrics before correlating with EDSS, rather than running the model with uncalibrated metrics and an addition protocol variable. Previously we found that the calibration method is generally superior to the protocol method under certain estimated conditions (i.e. subject variance, scanner variance, number of healthy controls at each site). The second goal of this notebook is to test a different calibration method based on transforming the simulated brain volume metrics into percentiles.
We simulate data the same way we did in the previous analyses:
%pylab inline
from utils import *
from IPython.html.widgets import interact, interactive, fixed
from IPython.html import widgets
from IPython.display import clear_output, display, HTML
import sklearn.linear_model as lm
import time
from mpld3 import enable_notebook, disable_notebook
from mpld3 import plugins
params = {'n_hc':50, #number of controls
'n_ms':[250,250], #number of MS patients at each site
'age_hc_range':(20,52), # Age range of the controls
'age_ms_range':[(20,45),(20,45)], # Age range of MS patients at each site
'd_dist':[.65,.65], # Distribution weighting of EDSS scores for each site
'alpha_a':-3, # Real coefficient of age
'alpha_d':[-16,-16], # Real coefficient of EDSS score
'b0':1644, #Brain Volume Intercept (in ccs)
'sig_b':[57,57], # Brain variation across subjects (in ccs) for each site
'sig_s':[16,16], # Variability in brain volume measurement across scanners for each site
'alpha':[0.7,1.3], # Scaling of each scanner
'gamma':[20,-20], # Offset of each scanner
'lo_mean':[2,2], # EDSS distribution parameters at both sites
'hi_mean':[6,6], #EDSS distribution parameters at both sites
'lo_sig': [1.1,1.1], #EDSS distribution parameters at both sites
'hi_sig': [0.8,0.8], #EDSS distribution parameters at both sites
}
sim = sim_data(**params)
Here we see the Bimodal-Normal distribution of EDSS at both sites.
The parameters across both sites are identical, so we expect the distribution to be very similar.
enable_notebook()
sns.set_style("white")
_,x=subplots(1)
ppl.remove_chartjunk(x,["top","right"])
x.set_title("EDSS Distribution",size=14)
sns.kdeplot(sim["dis_1"],shade=True,ax=x,label="Site 1");
sns.kdeplot(sim["dis_2"],shade=True,ax=x, label="Site 2");
x.set_xlim([0,10])
x.legend();
Next we will show the multivariate relationship between brain volume, age, and EDSS scores. The graphs show contours along which brain volume is constant. We expect the slopes of the contours to be similar but not identical across sites:
_,x = subplots(1,2,figsize=(12,4))
sns.set_style("darkgrid")
df1 = pd.DataFrame(index=range(params["n_ms"][0]))
df1["BV"] = sim["BV_ms_1"]
df1["age"] = sim["age_ms_1"]
df1["EDSS"] = sim["dis_1"]
df1["Group"] = np.array(["Lo","Hi"])[(sim["llr_1"] > 0).astype(int)]
iax = sns.interactplot("age","EDSS","BV", df1,cmap="coolwarm", filled=True, levels=25,ax=x[0],colorbar=True);
#iax.figure.colorbar(iax)
sns.set_style("darkgrid")
df2 = pd.DataFrame(index=range(params["n_ms"][1]))
df2["BV"] = sim["BV_ms_2"]
df2["age"] = sim["age_ms_2"]
df2["EDSS"] = sim["dis_2"]
df2["Group"] = np.array(["Lo","Hi"])[(sim["llr_2"] > 0).astype(int)]
sns.interactplot("age","EDSS","BV", df2,cmap="coolwarm", filled=True, levels=25,ax=x[1],colorbar=True);
x[1].set_title("Site 2");
x[0].set_title("Site 1");
For both sites, we see that for a given brain volume, as age increases, EDSS decreases. But the slopes do not match perfectly. This is because of the difference in scanner scaling factors.
The Log Likelihood Ratio (LLR) is defined as:
Where \(Pr_{\theta_L}\) is the probability density function(PDF) of a normal distribution centered at an EDSS 2 with a standard deviation of 1.1 (the low group), and \(Pr_{\theta_H}\) is the PDF of a normal distribution centered at EDSS of 6 with a standard deviation of 0.8.
We define the groups as:
Low: \(LLR > 0\)
High: \(LLR < 0\)
Below we plot the distribution of log likelihood ratios at each site:
_,x = subplots(1,2,figsize=(10,5))
[ppl.remove_chartjunk(i,["top","right"]) for i in x]
x[1].set_xlabel("EDSS LLR")
x[1].set_title("Site 2",size=14)
x[0].set_xlabel("EDSS LLR")
x[0].set_title("Site 1",size=14)
ppl.hist(x[0],sim["llr_1"],color=ppl.set1[0],alpha=0.7);
ppl.hist(x[1],sim["llr_2"],color=ppl.set1[1],alpha=0.7);
suptitle("EDSS Lo/Hi Log-Likelihood Ratio", size=16);
Next we will show the plots of brain volume, age and EDSS with the Low and High groups:
_,x=subplots(2,2,figsize=(12,12))
x = x.ravel()
llr = sim["llr_1"]
ppl.plot(x[0],sim["age_ms_1"][llr],sim["BV_ms_1"][llr],marker='o',linewidth=0,
markersize=10,color=ppl.set2[0],label="Lo",mew=0.5,mec=ppl.almost_black,alpha=0.9);
ppl.plot(x[0],sim["age_ms_1"][llr==False],sim["BV_ms_1"][llr==False],marker='o',linewidth=0,
markersize=10,color=ppl.set2[1],label="Hi",mew=0.5,mec=ppl.almost_black,alpha=0.9);
ppl.plot(x[1],sim["dis_1"][llr],sim["BV_ms_1"][llr],marker='o',linewidth=0,
markersize=10,color=ppl.set2[0],label="Lo",mew=0.5,mec=ppl.almost_black,alpha=0.9);
ppl.plot(x[1],sim["dis_1"][llr==False],sim["BV_ms_1"][llr==False],marker='o',linewidth=0,
markersize=10,color=ppl.set2[1],label="Hi",mew=0.5,mec=ppl.almost_black,alpha=0.9);
x[0].set_ylabel("Brain Volume",size=14)
x[0].set_xlabel("Age",size=14)
x[1].set_xlabel("EDSS",size=14)
llr = sim["llr_2"]
ppl.plot(x[2],sim["age_ms_2"][llr],sim["BV_ms_2"][llr],marker='o',linewidth=0,
markersize=10,color=ppl.set2[0],label="Lo",mew=0.5,mec=ppl.almost_black,alpha=0.9);
ppl.plot(x[2],sim["age_ms_2"][llr==False],sim["BV_ms_2"][llr==False],marker='o',linewidth=0,
markersize=10,color=ppl.set2[1],label="Hi",mew=0.5,mec=ppl.almost_black,alpha=0.9);
ppl.plot(x[3],sim["dis_2"][llr],sim["BV_ms_2"][llr],marker='o',linewidth=0,
markersize=10,color=ppl.set2[0],label="Lo",mew=0.5,mec=ppl.almost_black,alpha=0.9);
ppl.plot(x[3],sim["dis_2"][llr==False],sim["BV_ms_2"][llr==False],marker='o',linewidth=0,
markersize=10,color=ppl.set2[1],label="Hi",mew=0.5,mec=ppl.almost_black,alpha=0.9);
x[2].set_ylabel("Brain Volume",size=14)
x[2].set_xlabel("Age",size=14)
x[3].set_xlabel("EDSS",size=14)
ppl.legend(x[2]);
x[2].set_title("Site 2", size=14)
x[3].set_title("Site 2", size=14)
x[1].set_title("Site 1", size=14)
x[0].set_title("Site 1", size=14)
ppl.legend(x[3]);
ppl.legend(x[0]);
#suptitle("Site 1",size=14)
ppl.legend(x[1]);
suptitle("EDSS Severity Groups",size=16);
For each site, we convert the brain volume metric to percentile by calculating the empirical cumulative distribution function (CDF). The percentile is the CDF evaluated at a given brain volume.
from statsmodels.distributions.empirical_distribution import ECDF
def to_percentile(sim,ref_key,key):
foo = ECDF(sim[ref_key],side="right")
#get_p = lambda foo,sim,i : foo.y[np.nonzero(foo.x == sim[key][i])[0]][0]
#p = [get_p(foo,sim,i) for i in range(len(sim[key]))]
return foo(sim[key])
sim["pBV_ms_1"] = to_percentile(sim,"BV_hc_1","BV_ms_1")
sim["pBV_ms_2"] = to_percentile(sim,"BV_hc_2","BV_ms_2")
_,x=subplots(1,2,figsize=(12,4))
ppl.remove_chartjunk(x[0],["top","right"])
ppl.remove_chartjunk(x[1],["top","right"])
sns.kdeplot(to_percentile(sim,"BV_hc_1","BV_ms_1"),ax=x[0])
sns.kdeplot(to_percentile(sim,"BV_hc_2","BV_ms_2"),ax=x[0])
sns.kdeplot(sim["BV_ms_1"],ax=x[1])
sns.kdeplot(sim["BV_ms_2"],ax=x[1])
x[0].set_title("Percentile BV MS",size=14);
x[0].set_xlim([0,1])
x[1].set_title("Uncalibrated BV MS",size=14);
age_euro = np.array([24,24,24,38,37,32,28,30,32,32,50])
n_euro = len(age_euro)
BV_euro_real = params['alpha_a']*age_euro + np.random.normal(params['b0'],params['sig_b'][0],n_euro)
sns.set_style("white")
# Lets get scanned twice
#BV_euro_real = np.hstack((BV_euro_real,BV_euro_real))
#age_euro = np.hstack((age_euro,age_euro))
#n_euro = len(age_euro)
sim['BV_euro_1'] = params['alpha'][0]*BV_euro_real+np.random.normal(params['gamma'][0],params['sig_s'][0]*1/sqrt(2),n_euro)
sim['BV_euro_2'] = params['alpha'][1]*BV_euro_real+np.random.normal(params['gamma'][1],params['sig_s'][1]*1/sqrt(2),n_euro)
fig,ax = subplots(1)
sns.kdeplot(sim['BV_euro_1'],shade=True,label="Site 1 Euro",color=ppl.set1[0],ax=ax)
sns.kdeplot(sim['BV_euro_2'],shade=True,label="Site 2 Euro",color=ppl.set1[1],ax=ax)
xlabel("Brain Volume")
title("Brain Volume Distribution of EUROTRIP Subjects",size=14);
ppl.remove_chartjunk(ax,["top","right"])
sns.set_style("darkgrid")
fig,x = subplots(1,2,figsize=(12,4))
"""
ppl.plot(x[0],sim["BV_hc_1"],to_percentile(sim,"BV_hc_1","BV_hc_1"),marker='o',
linewidth=0,label="HC 1",markersize=15,alpha=0.5,mew=0.8,color=ppl.set2[0])
ppl.plot(x[0],sim["BV_ms_1"],to_percentile(sim,"BV_hc_1","BV_ms_1"),
marker='o',linewidth=0,label="MS 1",markersize=5,alpha=1,mew=0.5,color=ppl.almost_black)"""
ppl.scatter(x[0],sim["BV_hc_1"],to_percentile(sim,"BV_hc_1","BV_hc_1"),
marker='o',label="HC 1",alpha=0.5, color=ppl.set2[0],s=200)
ppl.scatter(x[0],sim["BV_ms_1"],to_percentile(sim,"BV_hc_1","BV_ms_1"),
marker='o',label="MS 1",alpha=1,color=ppl.almost_black,s=25)
Y1 = to_percentile(sim,"BV_hc_1","BV_euro_1")
#ppl.plot(x[0],sim["BV_euro_1"],Y1,
# marker='o',linewidth=0,label="msHC 1",markersize=10,alpha=1,mew=0.5,color=ppl.set2[1])
scatter1 = ppl.scatter(x[0],sim["BV_euro_1"],Y1,
marker='o',label="msHC 1",alpha=1,color=ppl.set2[1],s=100)
"""ppl.plot(x[1],sim["BV_hc_2"],to_percentile(sim,"BV_hc_2","BV_hc_2"),
marker='o',linewidth=0,label="HC 2",markersize=15,alpha=0.5,mew=0.5,color=ppl.set2[0])
ppl.plot(x[1],sim["BV_ms_2"],to_percentile(sim,"BV_hc_2","BV_ms_2"),
marker='o',linewidth=0,label="MS 2",markersize=5,alpha=1,mew=0.5,color=ppl.almost_black)"""
ppl.scatter(x[1],sim["BV_hc_2"],to_percentile(sim,"BV_hc_2","BV_hc_2"),
marker='o',label="HC 2",alpha=0.5, color=ppl.set2[0],s=200)
ppl.scatter(x[1],sim["BV_ms_2"],to_percentile(sim,"BV_hc_2","BV_ms_2"),
marker='o',label="MS 2",alpha=1,color=ppl.almost_black,s=25)
Y2 = to_percentile(sim,"BV_hc_2","BV_euro_2")
#ppl.plot(x[1],sim["BV_euro_2"],Y2,
# marker='o',linewidth=0,label="msHC 2",markersize=10,alpha=1,mew=0.5,color=ppl.set2[1])
scatter2 = ppl.scatter(x[1],sim["BV_euro_2"],Y2,
marker='o',label="msHC 2",alpha=1,color=ppl.set2[1],s=100)
img_str = '%d: %2.2f, %2.2f'
labels = [img_str%(idx,X,Y2[idx]) for idx,X in enumerate(sim["BV_euro_2"])]
tooltip = plugins.PointLabelTooltip(scatter2, labels)
plugins.connect(fig, tooltip);
labels = [img_str%(idx,X,Y1[idx]) for idx,X in enumerate(sim["BV_euro_1"])]
tooltip = plugins.PointLabelTooltip(scatter1, labels)
plugins.connect(fig, tooltip);
x[0].set_ylim([-0.1,1.1])
x[1].set_ylim([-0.1,1.1])
x[0].set_ylabel("Percentile",size=14)
x[0].set_xlabel("Brain Volume",size=14)
x[0].set_title("Site 1",size=16)
x[1].set_title("Site 2",size=16)
x[1].set_ylabel("Percentile",size=14)
x[1].set_xlabel("Brain Volume",size=14)
#ppl.legend(x[0],loc=0,fontsize=14);
#ppl.legend(x[1],loc=0,fontsize=14);
pdiff = np.abs(to_percentile(sim,"BV_hc_2","BV_euro_2")-to_percentile(sim,"BV_hc_1","BV_euro_1"))
p,t,c = test_all(pdiff,linspace(0.01,0.25,40))
fig,x=subplots(1,3,figsize=(12,4))
m = np.mean(pdiff)
s = np.std(pdiff)
sns.kdeplot(pdiff,ax=x[0])
x[0].vlines(m,0,5)
x[0].vlines(m+s,0,5,linestyle='--',color="grey")
x[0].vlines(m-s,0,5,linestyle='--',color="grey")
x[0].set_xlim([0,0.3]);
x[0].set_ylim([0,5]);
x[0].set_title("Percentile Differences for MultiSite Controls",size=14)
x[1].set_title("Power Curve",size=14)
x[2].set_title("Equivalence Test Cutoff",size=14)
ppl.plot(x[2],linspace(0.01,0.25,40),t,color=ppl.set2[0],linewidth=3,linestyle="--",label="T")
ppl.plot(x[2],linspace(0.01,0.25,40),c,color=ppl.set2[1],linewidth=3,label="Cutoff")
ppl.legend(x[2],fontsize=14,loc=0)
powercurve = ppl.scatter(x[1],linspace(0.01,0.25,40)[:20],p[:20],marker='o',s=100);
img_str = '%2.2f, %2.2f'
labels = [img_str%(X,p[:20][idx]) for idx,X in enumerate(linspace(0.01,0.25,40)[:20])]
tooltip = plugins.PointLabelTooltip(powercurve, labels)
plugins.connect(fig, tooltip);
from IPython.parallel.client.client import Client
from itertools import repeat
c = Client()
view = c[:]
view.block = True
def eq_test(y,ep):
import rpy2.robjects as ro
eq = ro.packages.importr("equivalence")
import numpy as np
x = ro.FloatVector(y)
e = ep/np.std(y)
res = eq.ptte_data(x,Epsilon=e)
P = res.rx("Power")[0][0]
T = res.rx("Tstat")[0][0]
C = res.rx("cutoff")[0][0]
return P,T,C
def test_all_parallel(pdiff,test_E = arange(0.001,0.01,0.0005),view=None):
g = view.map(eq_test,[pdiff]*len(test_E),test_E)
g = np.asarray(g)
return g
We first run a logistic regression of the two sites independently: with brain volume and age and try to predict group membership:
Site 1:
import sklearn.linear_model as lm
a = lm.LogisticRegression()
X1 = np.vstack((sim["BV_ms_1"],sim["age_ms_1"])).T
y1 = sim["llr_1"]>0
X1 = sm.add_constant(X1)
a.fit(X1,y1)
acc1 = a.score(X1,y1)
print acc1
Site 2:
a = lm.LogisticRegression()
X2 = np.vstack((sim["BV_ms_2"],sim["age_ms_2"])).T
X2 = sm.add_constant(X2)
y2 = sim["llr_2"]>0
a.fit(X2,sim["llr_2"]>0)
acc2 = a.score(X2,sim["llr_2"]>0)
print acc2
_,x = subplots(1,2,figsize=(12,4))
sns.set_style("darkgrid")
df1 = pd.DataFrame(index=range(params["n_ms"][0]))
df1["pBV"] = sim["pBV_ms_1"]
df1["age"] = sim["age_ms_1"]
df1["EDSS"] = sim["dis_1"]
df1["Low Group"] = sim["llr_1"] > 0
iax = sns.interactplot("age","pBV","Low Group", df1,
cmap="coolwarm", filled=True, levels=25,ax=x[0],
colorbar=True,logistic=True,
scatter_kws={"color": "black"},
contour_kws={"alpha": .8});
#iax.figure.colorbar(iax)
sns.set_style("darkgrid")
df2 = pd.DataFrame(index=range(params["n_ms"][1]))
df2["pBV"] = sim["pBV_ms_2"]
df2["age"] = sim["age_ms_2"]
df2["EDSS"] = sim["dis_2"]
df2["Low Group"] = sim["llr_2"] > 0
sns.interactplot("age","pBV","Low Group", df2,
cmap="coolwarm", filled=True, levels=25,ax=x[1],
colorbar=True,logistic=True,
scatter_kws={"color": "black"},
contour_kws={"alpha": .8});
x[1].set_title("Site 2\nAccuracy = %2.2f"%(acc2),size=16);
x[0].set_title("Site 1\n Accuracy = %2.2f"%acc1,size=16);
We then calculate and plot the ROC curve for both sites:
def get_ROC(X,y):
from sklearn.metrics import roc_curve, auc
from sklearn.utils import shuffle
import sklearn.linear_model as lm
X, y = shuffle(X, y)
half = int(len(y) / 2)
X_train, X_test = X[:half], X[half:]
y_train, y_test = y[:half], y[half:]
# Run classifier
classifier = lm.LogisticRegression()
classifier.fit(X_train, y_train)
probas_ = classifier.predict_proba(X_test)
# Compute ROC curve and area the curve
fpr, tpr, thresholds = roc_curve(y_test, probas_[:, 1])
roc_auc = auc(fpr, tpr)
return fpr,tpr,roc_auc
roc_info = view.map(get_ROC,[X1]*100,[y1]*100)
roc_info2 = view.map(get_ROC,[X2]*100,[y2]*100)
roc_info = np.asarray(roc_info)
roc_info2 = np.asarray(roc_info2)
fpr = roc_info[:,0]
tpr = roc_info[:,1]
auc = roc_info[:,2]
fpr2 = roc_info2[:,0]
tpr2 = roc_info2[:,1]
auc2 = roc_info2[:,2]
ROC Curve:
_,x=subplots(1,2,figsize=(12,4))
x[0].errorbar(np.mean(fpr),np.mean(tpr),xerr = np.std(fpr),yerr=np.std(tpr),alpha=0.15,color=ppl.set1[0],linewidth=3)
x[0].plot(np.mean(fpr),np.mean(tpr),color=ppl.set1[0],linewidth=3,label="Site 1")
x[0].errorbar(np.mean(fpr2),np.mean(tpr2),xerr = np.std(fpr2),yerr=np.std(tpr2),color=ppl.set1[1],alpha=0.15,linewidth=3)
x[0].plot(np.mean(fpr2),np.mean(tpr2),color=ppl.set1[1],linewidth=3,label="Site 2")
x[0].set_xlim([0,1])
x[0].set_ylim([0,1])
ppl.plot(x[0],[0, 1], [0, 1],color=ppl.set1[2],linewidth=2,linestyle='--')
x[1].set_title("Area Under the Curve",size=16)
x[0].set_title("ROC Curve",size=16)
x[0].set_xlabel("False Positive Rate")
x[0].set_xlabel("True Positive Rate")
x[0].legend(loc=0,fontsize=14)
ppl.boxplot(x[1],[auc,auc2]);
x[1].set_xticklabels(["Site 1","Site 2"]);
Next, we combine the data, run the regression and compute the ROC curve.
df_all = pd.DataFrame(index=range(sum(params["n_ms"])))
df_all["age"] = sim["age_all_ms"]
df_all["pBV"] = np.hstack((sim["pBV_ms_1"],sim["pBV_ms_2"]))
df_all["llr"] = np.hstack((sim["llr_1"],sim["llr_2"]))>0
a = lm.LogisticRegression()
X = np.vstack((df_all["pBV"],df_all["age"])).T
X = sm.add_constant(X)
y = df_all["llr"]
a.fit(X,y)
acc = a.score(X,y)
roc_info_all = view.map(get_ROC,[X]*100,[y]*100)
roc_info_all = np.asarray(roc_info_all)
fpr_all = roc_info_all[:,0]
tpr_all = roc_info_all[:,1]
auc_all = roc_info_all[:,2]
minlen = min([len(x) for x in fpr_all])
fpr_all = np.asarray([x[:minlen] for x in fpr_all])
tpr_all = np.asarray([x[:minlen] for x in tpr_all])
# Plot ROC curve
_,x=subplots(1,2,figsize=(12,4))
sns.interactplot("age","pBV","llr", df_all,cmap="coolwarm",
filled=True, levels=20,colorbar=True,logistic=True,ax=x[0],
scatter_kws={"color": "black"},
contour_kws={"alpha": .8});
x[0].set_title("Logistic Regression\n %2.2f Accuracy"%(acc),size= 16)
x[1].errorbar(np.mean(fpr),np.mean(tpr),xerr = np.std(fpr),yerr=np.std(tpr),alpha=0.15,color=ppl.set1[0],linewidth=3)
x[1].plot(np.mean(fpr),np.mean(tpr),color=ppl.set1[0],linewidth=3,label="Site 1, AUC = %2.2f"%(np.mean(auc)))
x[1].errorbar(np.mean(fpr2),np.mean(tpr2),xerr = np.std(fpr2),yerr=np.std(tpr2),color=ppl.set1[1],alpha=0.15,linewidth=3)
x[1].plot(np.mean(fpr2),np.mean(tpr2),color=ppl.set1[1],linewidth=3,label="Site 2, AUC = %2.2f"%(np.mean(auc2)))
x[1].errorbar(np.mean(fpr_all,axis=0),np.mean(tpr_all,axis=0),xerr = np.std(fpr_all,axis=0),
yerr=np.std(tpr_all,axis=0),alpha=0.15,color=ppl.set1[2],linewidth=3)
x[1].plot(np.mean(fpr_all,axis=0),np.mean(tpr_all,axis=0),color=ppl.set1[2],linewidth=3,
label="Combined, AUC = %2.2f"%(np.mean(auc_all)))
ppl.plot(x[1],[0, 1], [0, 1],color=ppl.set1[3],linewidth=2,linestyle="--")
x[1].set_xlim([0.0, 1.0])
x[1].set_ylim([0.0, 1.0])
x[1].set_xlabel('False Positive Rate', size=14)
x[1].set_ylabel('True Positive Rate',size=14)
x[1].set_title('ROC',size=16)
x[1].legend(loc="lower right",fontsize=14);
Here we see that the ROC curve has improved with combined, calibrated data
We use the multisite controls to test the percentile calibration. We need to find what change in percentile brain volume to use as an equivalence margin \(\delta\) to distinguish biased HC populations.
params = {'n_hc':50, #number of controls
'n_ms':[250,250], #number of MS patients at each site
'age_hc_range':(20,52), # Age range of the controls
'age_ms_range':[(20,45),(20,45)], # Age range of MS patients at each site
'd_dist':[.65,.65], # Distribution weighting of EDSS scores for each site
'alpha_a':-3, # Real coefficient of age
'alpha_d':[-16,-16], # Real coefficient of EDSS score
'b0':1644, #Brain Volume Intercept (in ccs)
'sig_b':[57,57], # Brain variation across subjects (in ccs) for each site
'sig_s':[16,16], # Variability in brain volume measurement across scanners for each site
'alpha':[0.7,1.3], # Scaling of each scanner
'gamma':[20,-20], # Offset of each scanner
'lo_mean':[2,2], # EDSS distribution parameters at both sites
'hi_mean':[6,6], #EDSS distribution parameters at both sites
'lo_sig': [1.1,1.1], #EDSS distribution parameters at both sites
'hi_sig': [0.8,0.8], #EDSS distribution parameters at both sites
}
b = boot(n_boot=100,keep_subs=False,params = params,age_euro_orig = age_euro,n_scans=1,sim_type=run_percentile_sim)
foo = [test_all_parallel(diff,linspace(0.01,0.2,20),view) for diff in b]
foo = np.asarray(foo)
_,x=subplots(1)
#pal = sns.blend_palette(["seagreen", "lightblue"])
sns.set_style("darkgrid")
pal = sns.color_palette("RdPu_r", 20)
sns.boxplot(foo[:,:,0],positions=range(20),ax=x,color=pal);
x.set_ylim([0,1.1])
x.set_title("Bootstrapped Power Analysis for Percentile Calibration",size=14)
x.set_xlabel("$\delta$ : Percentile",size=14)
x.set_ylabel("Power",size=14)
x.set_xticklabels(['%2.1f'%s for s in linspace(0.01,0.2,20)*100],rotation=90);
We only have the power to detect a 10% change in percentile across sites. This is really very high!
a = deepcopy(params)
n_boot = 100
age_euro = np.array([51,32,24,30,32,37,39,24,24,37,29,32,40,41])
alpha_enum = [[0.9,1.1],[0.8,1.2],[0.7,1.3]]
accept_diff = linspace(0.05,0.2,30)
Pall = np.zeros((len(alpha_enum),n_boot,len(accept_diff)))
Tall = np.zeros((len(alpha_enum),n_boot,len(accept_diff)))
Call = np.zeros((len(alpha_enum),n_boot,len(accept_diff)))
for i,sig in enumerate(alpha_enum):
a["alpha"] = sig
b = boot(n_boot=n_boot,keep_subs=True,params = a,age_euro_orig = age_euro,n_scans=1,sim_type=run_percentile_sim)
all_vals = np.asarray([test_all_parallel(diff,accept_diff,view) for diff in b])
Pall[i,:,:] = all_vals[:,:,0]
Tall[i,:,:] = all_vals[:,:,1]
Call[i,:,:] = all_vals[:,:,2]
average_pass = np.mean(Call>Tall,axis=1)
vector = [x[1] - x[0] for x in alpha_enum]
_,x = subplots(1)
for i,v in enumerate(vector):
ppl.plot(x,(np.array(accept_diff)+0.0*i)*100,average_pass[i],
label="alpha difference = %2.2f"%(v),color=ppl.set2[i],marker='o',markersize=10);
x.legend(loc=0)
x.set_xlabel("Similarity Margin(Difference in Percentile Brain Volume)",size=14)
x.set_ylabel("Pass Rate/ Similarity Probability",size=14)
x.set_ylim([0,1.1])
x.set_title("Probability of Similarity ",size=16)
#x.set_xlim([1,5]);
a0 = deepcopy(params)
n_boot = 100
age_euro = np.array([51,32,24,30,32,37,39,24,24,37,29,32,40,41])
n_hc = [25,50,75,100, 150,500]
accept_diff = linspace(0.05,0.2,30)
Pall = np.zeros((len(n_hc),n_boot,len(accept_diff)))
Tall = np.zeros((len(n_hc),n_boot,len(accept_diff)))
Call = np.zeros((len(n_hc),n_boot,len(accept_diff)))
t0 = time.time()
for i,hc in enumerate(n_hc):
a0["n_hc"] = hc
b = boot(n_boot=n_boot,keep_subs=True,params = a0,age_euro_orig = age_euro,n_scans=1,sim_type=run_percentile_sim)
#b = boot(n_boot=n_boot,keep_subs=True,params = a0,age_euro_orig = age_euro,sim_type=run_percentile_sim,n_scans=1)
all_vals = np.asarray([test_all_parallel(diff,accept_diff,view) for diff in b])
Pall[i,:,:] = all_vals[:,:,0]
Tall[i,:,:] = all_vals[:,:,1]
Call[i,:,:] = all_vals[:,:,2]
t_total = time.time()-t0
average_pass = np.mean(Call>Tall,axis=1)
vector = n_hc
_,x = subplots(1)
for i,v in enumerate(vector):
ppl.plot(x,(np.array(accept_diff)+0*i)*100,average_pass[i],
label="n hc = %d"%(v),color=ppl.set2[i],marker='o',markersize=10);
x.legend(loc=0)
x.set_xlabel("Similarity Margin(Difference in Percentile Brain Volume)",size=14)
x.set_ylabel("Pass Rate/ Similarity Probability",size=14)
x.set_ylim([0,1.1])
x.set_title("Probability of Similarity ",size=16)
#x.set_xlim([0.1,6]);
a = deepcopy(params)
n_boot = 100
age_euro = np.array([51,32,24,30,32,37,39,24,24,37,29,32,40,41])
n_scans = [1,2,3,4,5]
accept_diff = linspace(0.05,0.2,30)
Pall = np.zeros((len(n_scans),n_boot,len(accept_diff)))
Tall = np.zeros((len(n_scans),n_boot,len(accept_diff)))
Call = np.zeros((len(n_scans),n_boot,len(accept_diff)))
t0 = time.time()
for i,ns in enumerate(n_scans):
b = boot(n_boot=n_boot,keep_subs=True,params = a,age_euro_orig = age_euro,n_scans=ns, sim_type=run_percentile_sim)
all_vals = np.asarray([test_all_parallel(diff,accept_diff,view) for diff in b])
Pall[i,:,:] = all_vals[:,:,0]
Tall[i,:,:] = all_vals[:,:,1]
Call[i,:,:] = all_vals[:,:,2]
t_total = time.time()-t0
average_pass = np.mean(Call>Tall,axis=1)
vector = n_scans
_,x = subplots(1)
for i,v in enumerate(vector):
ppl.plot(x,(np.array(accept_diff)+0*i)*100,average_pass[i],
label="n scans = %d"%(v),color=ppl.set2[i],marker='o',markersize=10);
x.legend(loc=0)
x.set_xlabel("Similarity Margin( Percentile Difference in Brain Volume)",size=14)
x.set_ylabel("Pass Rate/ Similarity Probability",size=14)
x.set_ylim([0,1.1])
x.set_title("Probability of Similarity ",size=16)
#x.set_xlim([0.1,6]);
Bias_All = []
age_euro = np.array([51,32,24,30,32,37,39,24,24,37,29,32,40,41])
a_hc = [[1,1],[1,1.01],[1,1.025],[1,1.04],[1,1.05]]
for idx,a in enumerate(a_hc):
tmp = deepcopy(params)
tmp["alpha_hc"] = a
tmp["gamma_hc"] = [0, 0]
b = boot(n_boot=150,keep_subs=False,params = tmp,age_euro_orig = age_euro,n_scans=2,sim_type=run_percentile_sim)
foo1 = np.asarray([test_all_parallel(diff,linspace(0.05,0.4,50),view) for diff in b])
Bias_All.append(foo1)
_,x=subplots(1,figsize=(8,6))
b = [0,1,2.5,4,5]
for idx, B in enumerate(Bias_All):
ppl.plot(x,linspace(0.05,0.4,50)*100,np.mean(B[:,:,2]>B[:,:,1],axis=0),
marker='o',label="%2.1f Percent Bias"%(b[idx]),color=ppl.set2[idx],markersize=10)
#x.set_xticks(range(20))
x.set_xlabel("Equivalence Margin $\delta$ (\% BV)",size=14)
x.set_ylabel("Probability to detect Similarity",size=14)
x.legend(loc=0)
x.set_ylim([0,1.1]);
B_0 = np.mean(Bias_All[0][:,:,2]>Bias_All[0][:,:,1],axis=0)
B_1 = np.mean(Bias_All[1][:,:,2]>Bias_All[1][:,:,1],axis=0)
B_2 = np.mean(Bias_All[2][:,:,2]>Bias_All[2][:,:,1],axis=0)
#B_3 = np.mean(Bias_All[3][:,:,2]>Bias_All[3][:,:,1],axis=0)
plot(B_1,B_0,label="1 percent bias, AUC = %2.2f"%np.trapz(B_0,B_1),color=ppl.set2[1],linewidth=4);
plot(B_2,B_0,label="2.5 percent bias, AUC = %2.2f"%np.trapz(B_0,B_2),color=ppl.set2[2],linewidth=4);
#plot(B_3,B_0,label="4 percent bias, AUC = %2.2f"%np.trapz(B_0,B_3),color=ppl.set2[3],linewidth=4);
plot([0,1],[0,1],linestyle='--',color=ppl.almost_black,alpha=0.5);
title("ROC Curve",size=16)
xlabel("False Positive Rate",size=14)
ylabel("True Positive Rate",size=14)
ylim([-0.05,1.05])
xlim([-0.05,1.05])
legend(loc=4,fontsize=14);
The percentile calibration method is more sensitive to bias than the OLS method, as shown by the area under the ROC curves for given levels of bias
With the percentile calibration, we still see that we need >50 HCs and >2 scans per multisite control
We have the power to detect a 10 percentile equivalence margin, but the data is only equivalent at margins > 15 percentiles